/*
 * com.javable.cese -- Cell Electrophysiology Simulation Environment (CESE)
 *
 * Copyright (C) 1999 - 2005 Sergey Missan
 *
 * Dalhousie Univ., Halifax, NS, Canada
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */

package models.com.javable.cese.models;

import java.io.Serializable;

import com.javable.cese.templates.AbstractAPModel;

/**
 * The Faber-Rudy I Model of Mammalian Ventricular Cardiac Action
 * Potentials, 2007. Detailed list of equations and model description are
 * provided in
 * <p>
 * Circ Res 1991;68:1501-1526 <br>
 * Circ Res 1994;74:1071-1096 <br>
 * Circ Res 1994;74:1097-1113 <br>
 * Circ Res 1995;77:140-152 <br>
 * Biophys J 1995;68:949-964 <br>
 * Cardiovasc Res 1997;35:256-272 <br>
 * Circulation 1999;99:2466-2474 <br>
 * Cardiovas Res 1999;42:530-542 <br>
 * Nature 1999;400:566-569 <br>
 * Circulation 2000;101:1192-1198 <br>
 * Biophys J 2000;78:2392-2404 <br>
 * Biophys J 2007;92:1522-1543.
 */
public class FaberRudyModel_I_2007 extends AbstractAPModel implements
        Serializable {

    /* The beginning here lists all the variables and defines them */

    /* Cell Geometry */
    private double l = 0.01; // Length of the cell (cm)

    public synchronized void setl(double sl) {
        l = sl;
        getCellGeometry();
    }

    public double getl() {
        return l;
    }

    private double a = 0.0011; // Radius of the cell (cm)

    public synchronized void seta(double sa) {
        a = sa;
        getCellGeometry();
    }

    public double geta() {
        return a;
    }

    private double vcell; // Cell volume (uL)

    public double getvcell() {
        return vcell;
    }

    private double ageo; // Geometric membrane area (cm^2)

    public double getageo() {
        return ageo;
    }

    private double acap; // Capacitive membrane area (cm^2)

    public double getacap() {
        return acap;
    }

    private double vmyo; // Myoplasm volume (uL)

    public double getvmyo() {
        return vmyo;
    }

    private double vmito; // Mitochondria volume (uL)

    public double getvmito() {
        return vmito;
    }

    private double vsr; // SR volume (uL)

    public double getvsr() {
        return vsr;
    }

    private double vnsr; // NSR volume (uL)

    public double getvnsr() {
        return vnsr;
    }

    private double vjsr; // JSR volume (uL)

    public double getvjsr() {
        return vjsr;
    }

    //private double vcleft; // Cleft volume (uL)

    //public double getvcleft() {
    //    return vcleft;
    //}

    private double vss; // Subspace volume (uL)

    public double getvss() {
        return vss;
    }

    /* Terms for Solution of Conductance and Reversal Potential */
    public static final double R = 8314; // Universal Gas Constant (J/kmol*K)
    public static final double frdy = 96485; // Faraday's Constant (C/mol)

    private double temp = 310; // Temperature (K)

    public synchronized void settemp(double stemp) {
        temp = stemp;
    }

    public double gettemp() {
        return temp;
    }

    /* Ion Valences */
    public static final double zna = 1; // Na valence
    public static final double zk = 1; // K valence
    public static final double zca = 2; // Ca valence

    /* Ion Concentrations */
    private double nai = 10.71; // Intracellular Na Concentration (mM)

    public synchronized void setnai(double snai) {
        nai = snai;
    }

    public double getnai() {
        return nai;
    }

    private double nao = 130; // Extracellular Na Concentration (mM)

    public synchronized void setnao(double snao) {
        nao = snao;
    }

    public double getnao() {
        return nao;
    }

    private double nass = 10.71; // Subspace Na Concentration (mM)

    public synchronized void setnass(double snass) {
        nass = snass;
    }

    public double getnass() {
        return nass;
    }

    private double ki = 137.6; // Intracellular K Concentration (mM)

    public synchronized void setki(double ski) {
        ki = ski;
    }

    public double getki() {
        return ki;
    }

    private double ko = 4.5; // Extracellular K Concentration (mM)

    public synchronized void setko(double sko) {
        ko = sko;
    }

    public double getko() {
        return ko;
    }

    private double kss = 137.6; // Subspace K Concentration (mM)

    public synchronized void setkss(double skss) {
        kss = skss;
    }

    public double getkss() {
        return kss;
    }

    private double cai = 0.0000821; // Intracellular Ca Concentration (mM)

    public synchronized void setcai(double scai) {
        cai = scai;
    }

    public double getcai() {
        return cai;
    }

    private double cao = 1.8; // Extracellular Ca Concentration (mM)

    public synchronized void setcao(double scao) {
        cao = scao;
    }

    public double getcao() {
        return cao;
    }

    private double cass = 0.0000672; // Subspace Ca Concentration (mM)

    public synchronized void setcass(double scass) {
        cass = scass;
    }

    public double getcass() {
        return cass;
    }

    private double cmdn = 0.0024; // Calmodulin Buffered Ca Concentration (mM)

    public synchronized void setcmdn(double scmdn) {
        cmdn = scmdn;
    }

    public double getcmdn() {
        return cmdn;
    }

    private double trpn = 0.0133; // Troponin Buffered Ca Concentration (mM)

    public synchronized void settrpn(double strpn) {
        trpn = strpn;
    }

    public double gettrpn() {
        return trpn;
    }

    private double nsr = 2.553; // NSR Ca Concentration (mM) (1.8 mM in LRdII)

    public synchronized void setnsr(double snsr) {
        nsr = snsr;
    }

    public double getnsr() {
        return nsr;
    }

    private double jsr = 2.479; // JSR Ca Concentration (mM) (1.8 mM in LRdII)

    public synchronized void setjsr(double sjsr) {
        jsr = sjsr;
    }

    public double getjsr() {
        return jsr;
    }

    private double csqn = 6.9; // Calsequestrin Buffered Ca Concentration (mM)

    public synchronized void setcsqn(double scsqn) {
        csqn = scsqn;
    }

    public double getcsqn() {
        return csqn;
    }

    /* Myoplasmic Na Ion Concentration Changes */
    private double naiont; // Total Na Ion Flow (mM/ms)

    public double getnaiont() {
        return naiont;
    }

    private double dnai; // Change in Intracellular Na Concentration (mM)

    public double getdnai() {
        return dnai;
    }

    /* Myoplasmic K Ion Concentration Changes */
    private double kiont; // Total K Ion Flow (mM/ms)

    public double getkiont() {
        return kiont;
    }

    private double dki; // Change in Intracellular K Concentration (mM)

    public double getdki() {
        return dki;
    }

    /* NSR Ca Ion Concentration Changes */
    private double dnsr; // Change in [Ca] in the NSR (mM)

    public double getdnsr() {
        return dnsr;
    }

    private double iup; // Ca uptake from myo. to NSR (mM/ms)

    public double getiup() {
        return iup;
    }

    private double ileak; // Ca leakage from NSR to myo. (mM/ms)

    public double getileak() {
        return ileak;
    }

    private double kleak; // Rate constant of Ca leakage from NSR (ms^-1)

    public double getkleak() {
        return kleak;
    }

    private double kmup = 0.00092; // Half-saturation concentration of iup (mM)

    public synchronized void setkmup(double skmup) {
        kmup = skmup;
    }

    public double getkmup() {
        return kmup;
    }

    private double iupbar = 0.017325; // Max. current through iup channel (mM/ms), 0.005 in LRdII

    public synchronized void setiupbar(double siupbar) {
        iupbar = siupbar;
    }

    public double getiupbar() {
        return iupbar;
    }

    private double nsrbar = 15; // Max. [Ca] in NSR (mM)

    public synchronized void setnsrbar(double snsrbar) {
        nsrbar = snsrbar;
    }

    public double getnsrbar() {
        return nsrbar;
    }

    /* JSR Ca Ion Concentration Changes */
    private double djsr; // Change in [Ca] in the JSR (mM)

    public double getdjsr() {
        return djsr;
    }

    //private double tauon = 2;

    // Time constant of activation of Ca release from JSR (ms)
    //public synchronized void settauon(double stauon) {
    //    tauon = stauon;
    //}

    //public double gettauon() {
    //    return tauon;
    //}

    //private double tauoff = 2;

    // Time constant of deactivation of Ca release from JSR (ms)
    //public synchronized void settauoff(double stauoff) {
    //    tauoff = stauoff;
    //}

    //public double gettauoff() {
    //    return tauoff;
    //}

    //private double caitzero;

    // Total myoplasmic [Ca] at time of dv/dt max. (mM)
    //public double getcaitzero() {
    //    return caitzero;
    //}

    //private double caiatwo;

    // Total myoplasmic [Ca] 2 ms. after dv/dt max. (mM)
    //public double getcaiatwo() {
    //    return caiatwo;
    //}

    //private double dcaitwo;

    // Cumulative Ca entry into the cell 2 ms. after dv/dt max. (mM)
    //public double getdcaitwo() {
    //    return dcaitwo;
    //}

    //private double tcount = 3;
    // Counter for 2 ms. after onset of dv/dt max. (ms)

    //private double dcaith = 0.000242;

    // Threshold for external triggering of Ca release from JSR (mM)
    //public synchronized void setdcaith(double sdcaith) {
    //    dcaith = sdcaith;
    //}

    //public double getdcaith() {
    //    return dcaith;
    //}

    //private double grelbarcicr = 0;

    // Max. rate constant of Ca release from JSR due to CICR (ms^-1)
    //public synchronized void setgrelbarcicr(double sgrelbarcicr) {
    //    grelbarcicr = sgrelbarcicr;
    //}

    //public double getgrelbarcicr() {
    //    return grelbarcicr;
    //}

    //private double grelcicr;

    // Rate constant of Ca release from JSR due to CICR (ms^-1)
    //public double getgrelcicr() {
    //    return grelcicr;
    //}

    //private double tcicr; // t=0 at time of CICR (ms)

    //public double gettcicr() {
    //    return tcicr;
    //}

    //private double kmrel = 0.0008;

    // Half-saturation constant of irel for CICR (mM)
    //public synchronized void setkmrel(double skmrel) {
    //    kmrel = skmrel;
    //}

    //public double getkmrel() {
    //    return kmrel;
    //}

    //private double irelcicr; // Ca release from JSR to myo. due to CICR (mM/ms)

    //public double getirelcicr() {
    //    return irelcicr;
    //}

    //private double csqnth = 8.75;

    // Threshold for release of Ca from Calsequestrin (CSQN) due to JSR
    // overload (mM)
    //public synchronized void setcsqnth(double scsqnth) {
    //    csqnth = scsqnth;
    //}

    //public double getcsqnth() {
    //    return csqnth;
    //}

    //private double gmaxrel = 60;

    // Max. rate constant of Ca release from JSR due to overload (ms^-1)
    //public synchronized void setgmaxrel(double sgmaxrel) {
    //    gmaxrel = sgmaxrel;
    //}

    //public double getgmaxrel() {
    //    return gmaxrel;
    //}

    //private double grelbarjsrol = 0;

    // Rate constant of Ca release from JSR due to overload (ms^-1)
    //public synchronized void setgrelbarjsrol(double sgrelbarjsrol) {
    //    grelbarjsrol = sgrelbarjsrol;
    //}

    //public double getgrelbarjsrol() {
    //    return grelbarjsrol;
    //}

    //private double greljsrol;

    // Rate constant of Ca release from JSR due to JSR overload (ms^-1)
    //public double getgreljsrol() {
    //    return greljsrol;
    //}

    //private double tjsrol; // t=0 at time of JSR overload (ms)

    //public double gettjsrol() {
    //    return tjsrol;
    //}

    //private double ireljsrol;

    // Ca release from JSR to myo. due to JSR overload (mM/ms)
    //public double getireljsrol() {
    //    return ireljsrol;
    //}

    private double csqnbar = 10; // Max. [Ca] buffered in CSQN (mM)

    public synchronized void setcsqnbar(double scsqnbar) {
        csqnbar = scsqnbar;
    }

    public double getcsqnbar() {
        return csqnbar;
    }

    private double kmcsqn = 0.8;

    // Equilibrium constant of buffering for CSQN (mM)
    public synchronized void setkmcsqn(double skmcsqn) {
        kmcsqn = skmcsqn;
    }

    public double getkmcsqn() {
        return kmcsqn;
    }

    private double bjsr;

    // b Variable for analytical computation of [Ca] in JSR (mM)
    public double getbjsr() {
        return bjsr;
    }

    //private double cjsr;

    // c Variable for analytical computation of [Ca] in JSR (mM)
    //public double getcjsr() {
    //    return cjsr;
    //}

    /* Translocation of Ca Ions from NSR to JSR */
    private double itr;

    // Translocation current of Ca ions from NSR to JSR (mM/ms)
    public double getitr() {
        return itr;
    }

    private double tautr = 120;//180 ms in LRdII

    // Time constant of Ca transfer from NSR to JSR (ms)
    public synchronized void settautr(double stautr) {
        tautr = stautr;
    }

    public double gettautr() {
        return tautr;
    }

    /* Myoplasmic Ca Ion Concentration Changes */
    private double caiont; // Total Ca Ion Flow (mM/ms)

    public double getcaiont() {
        return caiont;
    }

    private double dcai; // Change in myoplasmic Ca concentration (mM)

    public double getdcai() {
        return dcai;
    }

    private double catotal; // Total myoplasmic Ca concentration (mM)

    public double getcatotal() {
        return catotal;
    }

    private double bmyo;

    // b Variable for analytical computation of [Ca] in myoplasm (mM)
    public double getbmyo() {
        return bmyo;
    }

    //private double cmyo;

    // c Variable for analytical computation of [Ca] in myoplasm (mM)
    //public double getcmyo() {
    //    return cmyo;
    //}

    //private double dmyo;

    // d Variable for analytical computation of [Ca] in myoplasm (mM)
    //public double getdmyo() {
    //    return dmyo;
    //}

    private double cmdnbar = 0.050;

    // Max. [Ca] buffered in Calmodulin (CMDN) (mM)
    public synchronized void setcmdnbar(double scmdnbar) {
        cmdnbar = scmdnbar;
    }

    public double getcmdnbar() {
        return cmdnbar;
    }

    private double trpnbar = 0.070;

    // Max. [Ca] buffered in Troponin (TRPN) (mM)catotal = trpn + cmdn + dcai + cai;
    public synchronized void settrpnbar(double strpnbar) {
        trpnbar = strpnbar;
    }

    public double gettrpnbar() {
        return trpnbar;
    }

    private double kmcmdn = 0.00238;

    // Equilibrium constant of buffering for CMDN (mM)
    public synchronized void setkmcmdn(double skmcmdn) {
        kmcmdn = skmcmdn;
    }

    public double getkmcmdn() {
        return kmcmdn;
    }

    private double kmtrpn = 0.0005;

    // Equilibrium constant of buffering for TRPN (mM)
    public synchronized void setkmtrpn(double skmtrpn) {
        kmtrpn = skmtrpn;
    }

    public double getkmtrpn() {
        return kmtrpn;
    }

    private double idiffca; //Ca transfer from subspace to myoplasm (mM/ms)

    public double getidiffca() {
        return idiffca;
    }

    private double bcass; // b Variable for analytical computation of [Ca] in subspace (mM)

    public double getbcass() {
        return bcass;
    }

    private double dcass; // Change in subspace Ca concentration (mM)

    public double getdcass() {
        return dcass;
    }

    private double bsrbar = 0.047;

    // Anionic SR Ca binding sites in subspace (mM)
    public synchronized void setbsrbar(double sbsrbar) {
        bsrbar = sbsrbar;
    }

    public double getbsrbar() {
        return bsrbar;
    }

    private double kmbsr = 0.00087;

    // Equilibrium constant for anionic SR Ca binding sites in subspace (mM)
    public synchronized void setkmbsr(double skmbsr) {
        kmbsr = skmbsr;
    }

    public double getkmbsr() {
        return kmbsr;
    }

    private double bslbar = 2.124;

    // Anionic sarcolemmal Ca binding sites in subspace (mM)
    public synchronized void setbslbar(double sbslbar) {
        bslbar = sbslbar;
    }

    public double getbslbar() {
        return bslbar;
    }

    private double kmbsl = 0.127;

    // Equilibrium constant for anionic sarcolemmal Ca binding sites in subspace (mM)
    public synchronized void setkmbsl(double skmbsl) {
        kmbsl = skmbsl;
    }

    public double getkmbsl() {
        return kmbsl;
    }

    private double taudiff = 0.1;

    // Time constant of Ca transfer from subspace to myoplasm (ms)
    public synchronized void settaudiff(double staudiff) {
        taudiff = staudiff;
    }

    public double gettaudiff() {
        return taudiff;
    }

    private double idiffna; //Na transfer from subspace to myoplasm (mM/ms)

    public double getidiffna() {
        return idiffna;
    }

    private double dnass; // Change in subspace Na concentration (mM)

    public double getdnass() {
        return dnass;
    }

    private double idiffk; //K transfer from subspace to myoplasm (mM/ms)

    public double getidiffk() {
        return idiffk;
    }

    private double dkss; // Change in subspace K concentration (mM)

    public double getdkss() {
        return dkss;
    }

    //private double gpig;
    // Tribute to all the guinea pigs killed for the advancement of knowledge

    /* Fast Sodium Current (time dependent) */
    private double ina; // Fast Na Current (uA/uF)

    public synchronized void setina(double sina) {
        ina = sina;
    }

    public double getina() {
        return ina;
    }

    private double gna; // Max. Conductance of the Na Channel (mS/uF)

    public synchronized void setgna(double sgna) {
        gna = sgna;
    }

    public double getgna() {
        return gna;
    }

    private double ena; // Reversal Potential of Na (mV)

    public synchronized void setena(double sena) {
        ena = sena;
    }

    public double getena() {
        return ena;
    }

    private double ah; // Na alpha-h rate constant (ms^-1)

    public double getah() {
        return ah;
    }

    private double bh; // Na beta-h rate constant (ms^-1)

    public double getbh() {
        return bh;
    }

    private double aj; // Na alpha-j rate constant (ms^-1)

    public double getaj() {
        return aj;
    }

    private double bj; // Na beta-j rate constant (ms^-1)

    public double getbj() {
        return bj;
    }

    private double am; // Na alpha-m rate constant (ms^-1)

    public double getam() {
        return am;
    }

    private double bm; // Na beta-m rate constant (ms^-1)

    public double gebm() {
        return bm;
    }

    private double h; // Na inactivation h

    public double geth() {
        return h;
    }

    private double j; // Na inactivation j

    public double getj() {
        return j;
    }

    private double m; // Na activation m

    public double getm() {
        return m;
    }

    private double scgna = 1;

    // Scaling factor for Max. Conductance of the Na Channel
    public synchronized void setscgna(double sscgna) {
        scgna = sscgna;
    }

    public double getscgna() {
        return scgna;
    }

    /* Current through L-type Ca Channel */
    private double ilca; // Ca current through L-type Ca channel (uA/uF)

    public synchronized void setilca(double silca) {
        ilca = silca;
    }

    public double getilca() {
        return ilca;
    }

    private double ilcana; // Na current through L-type Ca channel (uA/uF)

    public synchronized void setilcana(double silcana) {
        ilcana = silcana;
    }

    public double getilcana() {
        return ilcana;
    }

    private double ilcak; // K current through L-type Ca channel (uA/uF)

    public synchronized void setilcak(double silcak) {
        ilcak = silcak;
    }

    public double getilcak() {
        return ilcak;
    }

    private double ilcatot;

    // Total current through the L-type Ca channel (uA/uF)
    public double getilcatot() {
        return ilcatot;
    }

    private double ibarca; // Max. Ca current through Ca channel (uA/uF)

    public double getibarca() {
        return ibarca;
    }

    private double ibarna; // Max. Na current through Ca channel (uA/uF)

    public double getibarna() {
        return ibarna;
    }

    private double ibark; // Max. K current through Ca channel (uA/uF)

    public double getibark() {
        return ibark;
    }

    //newly added transition rates and states values for complex Markov model of ICaL (Faber&Rudy2007)
    private double icala; //ICaL activation rate alpha (ms^-1)

    public double geticala() {
        return icala;
    }

    private double icalb; //ICaL deactivation rate beta (ms^-1)

    public double geticalb() {
        return icalb;
    }

    private double icalgf; //ICaL fast VD inactivation rate form closed state gamma fast (ms^-1)

    public double geticalgf() {
        return icalgf;
    }

    private double icalgs; //ICaL slow VD inactivation rate form closed state gamma slow (ms^-1)

    public double geticalgs() {
        return icalgs;
    }

    private double icalff; //ICaL fast VD inactivation rate form open state fi fast (ms^-1)

    public double geticalff() {
        return icalff;
    }

    private double icalfs; //ICaL slow VD inactivation rate form open state fi slow (ms^-1)

    public double geticalfs() {
        return icalfs;
    }

    private double icallf; //ICaL fast VD recovery from inactivation rate in open state lambda fast (ms^-1)

    public double geticallf() {
        return icallf;
    }

    private double icalls; //ICaL slow VD recovery from inactivation rate in open state lambda slow (ms^-1)

    public double geticalls() {
        return icalls;
    }

    private double icalof; //ICaL fast VD recovery from inactivation rate in closed state omega fast (ms^-1)

    public double geticalof() {
        return icalof;
    }

    private double icalos; //ICaL slow VD recovery from inactivation rate in closed state omega slow (ms^-1)

    public double geticalos() {
        return icalos;
    }

    private double icalosf; //ICaL  switch rate from slow to fast inactivation  omega slow fast (ms^-1)

    public double geticalosf() {
        return icalosf;
    }

    private double icalofs; //ICaL  switch rate from fast to slow inactivation  omega fast slow (ms^-1)

    public double geticalofs() {
        return icalofs;
    }

    private double icald; //ICaL  switch rate from mode V to mode Ca delta (ms^-1)

    public double geticald() {
        return icald;
    }

    private double icalt = 0.01; 
    
    //ICaL  switch rate from mode Ca to mode V theta (ms^-1)
    public synchronized void seticalt(double scalt) {
        icalt = scalt;
    }

    private double icalc0 = 0.923; 
    
    //ICaL  probability of state C0
    public synchronized void seticalc0(double scalc0) {
        icalc0 = scalc0;
    }

    private double icalc1 = 0.049; 
    
    //ICaL  probability of state C1
    public synchronized void seticalc1(double scalc1) {
        icalc1 = scalc1;
    }

    private double icalc2 = 0; 
    
    //ICaL  probability of state C2
    public synchronized void seticalc2(double scalc2) {
        icalc2 = scalc2;
    }

    private double icalc3 = 0; 
    
    //ICaL  probability of state C3
    public synchronized void seticalc3(double scalc3) {
        icalc3 = scalc3;
    }

    private double icalo = 0; 
    
    //ICaL  probability of state O
    public synchronized void seticalo(double scalo) {
        icalo = scalo;
    }

    private double icalif = 0; 
    
    //ICaL  probability of state If
    public synchronized void seticalif(double scalif) {
        icalif = scalif;
    }

    private double icalis = 0; 
    
    //ICaL  probability of state Is
    public synchronized void seticalis(double scalis) {
        icalis = scalis;
    }

    private double icalc0ca = 0.025; 
    
    //ICaL  probability of state C0Ca
    public synchronized void seticalc0ca(double scalc0ca) {
        icalc0ca = scalc0ca;
    }

    private double icalc1ca = 0.003; 
    
    //ICaL  probability of state C1Ca
    public synchronized void seticalc1ca(double scalc1ca) {
        icalc1ca = scalc1ca;
    }

    private double icalc2ca = 0; 
    
    //ICaL  probability of state C2Ca
    public synchronized void seticalc2ca(double scalc2ca) {
        icalc2ca = scalc2ca;
    }

    private double icalc3ca = 0; 
    
    //ICaL  probability of state C3Ca
    public synchronized void seticalc3ca(double scalc3ca) {
        icalc3ca = scalc3ca;
    }

    private double icalica = 0; 
    
    //ICaL  probability of state ICa
    public synchronized void seticalica(double scalica) {
        icalica = scalica;
    }

    private double icalifca = 0; 
    
    //ICaL  probability of state IfCa
    public synchronized void seticalifca(double scalifca) {
        icalifca = scalifca;
    }

    private double icalisca = 0; 
    
    //ICaL  probability of state IsCa
    public synchronized void seticalisca(double scalisca) {
        icalisca = scalisca;
    }

    //old ICaL model parameters commented (LuoRudyII2000)
    //private double ad; // Ca alpha-d rate constant (ms^-1)

    //public double getad() {
    //    return ad;
    //}

    //private double bd; // Ca beta-d rate constant (ms^-1)

    //public double getbd() {
    //    return bd;
    //}

    //private double d; // Voltage dependent activation gate d

    //public double getd() {
    //    return d;
    //}

    //private double dss; // Steady-state value of activation gate d

    //public double getdss() {
    //    return dss;
    //}

    //private double taud; // Time constant of gate d (ms^-1)

    //public double gettaud() {
    //    return taud;
    //}

    //private double af; // Ca alpha-f rate constant (ms^-1)

    //public double getaf() {
    //    return af;
    //}

    //private double bf; // Ca beta-f rate constant (ms^-1)

    //public double getbf() {
    //    return bf;
    //}

    //private double f; // Voltage dependent inactivation gate f

    //public double getf() {
    //    return f;
    //}

    //private double fss; // Steady-state value of inactivation gate f

    //public double getfss() {
    //    return fss;
    //}

    //private double tauf; // Time constant of gate f (ms^-1)

    //public double gettauf() {
    //    return tauf;
    //}

    //private double fca; // Ca dependent inactivation gate

    //public double getfca() {
    //    return fca;
    //}

    private double kmca = 0.0006;

    // Half-saturation concentration of Ca channel (mM)
    public synchronized void setkmca(double skmca) {
        kmca = skmca;
    }

    public double getkmca() {
        return kmca;
    }

    private double pca = 0.001215; // Permeability of membrane to Ca (cm/s), 0.00054 in LRdII

    public synchronized void setpca(double spca) {
        pca = spca;
    }

    public double getpca() {
        return pca;
    }

    private double gacai = 0.01; // Activity coefficient of Ca in

    public synchronized void setgacai(double sgacai) {
        gacai = sgacai;
    }

    public double getgacai() {
        return gacai;
    }

    private double gacao = 0.341; // Activity coefficient of Ca out

    public synchronized void setgacao(double sgacao) {
        gacao = sgacao;
    }

    public double getgacao() {
        return gacao;
    }

    private double pna = 0.000001518; // Permeability of membrane to Na (cm/s), 0.000000675 in LRdII

    public synchronized void setpna(double spna) {
        pna = spna;
    }

    public double getpna() {
        return pna;
    }

    private double ganai = 0.75; // Activity coefficient of Na in

    public synchronized void setganai(double sganai) {
        ganai = sganai;
    }

    public double getganai() {
        return ganai;
    }

    private double ganao = 0.75; // Activity coefficient of Na out

    public synchronized void setganao(double sganao) {
        ganao = sganao;
    }

    public double getganao() {
        return ganao;
    }

    private double pk = 0.000000434; // Permeability of membrane to K (cm/s), 0.000000193 in LRdII

    public synchronized void setpk(double spk) {
        pk = spk;
    }

    public double getpk() {
        return pk;
    }

    private double gaki = 0.75; // Activity coefficient of K in

    public synchronized void setgaki(double sgaki) {
        gaki = sgaki;
    }

    public double getgaki() {
        return gaki;
    }

    private double gako = 0.75; // Activity coefficient of K out

    public synchronized void setgako(double sgako) {
        gako = sgako;
    }

    public double getgako() {
        return gako;
    }

    private double scgcal = 1;

    // Scaling factor for Max. Conductance of the L-type Ca channel
    public synchronized void setscgcal(double sscgcal) {
        scgcal = sscgcal;
    }

    public double getscgcal() {
        return scgcal;
    }

    /* Current through RyR Channel */
    private double irel; // Ca current through RyR channel (uA/uF)

    public synchronized void setirel(double sirel) {
        irel = sirel;
    }

    public double getirel() {
        return irel;
    }

    private double iryra1; //RyR activation rate alpha1 (ms^-1)

    public double getiryra1() {
        return iryra1;
    }

    private double iryra2; //RyR activation rate alpha2 (ms^-1)

    public double getiryra2() {
        return iryra2;
    }

    private double iryrb1; //RyR activation rate beta1 (ms^-1)

    public double getiryrb1() {
        return iryrb1;
    }

    private double iryrb2; //RyR activation rate beta2 (ms^-1)

    public double getiryrb2() {
        return iryrb2;
    }

    private double iryrb3; //RyR activation rate beta3 (ms^-1)

    public double getiryrb3() {
        return iryrb3;
    }

    private double iryrb4; //RyR activation rate beta4 (ms^-1)

    public double getiryrb4() {
        return iryrb4;
    }

    private double iryrg1; //RyR activation rate gamma1 (ms^-1)

    public double getiryrg1() {
        return iryrg1;
    }

    private double iryrg2; //RyR activation rate gamma2 (ms^-1)

    public double getiryrg2() {
        return iryrg2;
    }

    private double iryrg3; //RyR activation rate gamma3 (ms^-1)

    public double getiryrg3() {
        return iryrg3;
    }

    private double iryrg4; //RyR activation rate gamma4 (ms^-1)

    public double getiryrg4() {
        return iryrg4;
    }

    private double iryrg5; //RyR activation rate gamma5 (ms^-1)

    public double getiryrg5() {
        return iryrg5;
    }

    private double iryrd1; //RyR activation rate delta1 (ms^-1)

    public double getiryrd1() {
        return iryrd1;
    }

    private double iryrd2; //RyR activation rate delta2 (ms^-1)

    public double getiryrd2() {
        return iryrd2;
    }

    private double iryrd3; //RyR activation rate delta3 (ms^-1)

    public double getiryrd3() {
        return iryrd3;
    }

    private double iryrd4; //RyR activation rate delta4 (ms^-1)

    public double getiryrd4() {
        return iryrd4;
    }

    private double iryrd5; //RyR activation rate delta5 (ms^-1)

    public double getiryrd5() {
        return iryrd5;
    }

    private double iryrc1 = 0.89;

    //RyR  probability of state C1
    public synchronized void setiryrc1(double sryrc1) {
        iryrc1 = sryrc1;
    }

    private double iryrc2 = 0.021; 
    
    //RyR  probability of state C2
    public synchronized void setiryrc2(double sryrc2) {
        iryrc2 = sryrc2;
    }

    private double iryrc3 = 0.003; 
    
    //RyR  probability of state C3
    public synchronized void setiryrc3(double sryrc3) {
        iryrc3 = sryrc3;
    }

    private double iryrc4 = 0.001; 
    
    //RyR  probability of state C4
    public synchronized void setiryrc4(double sryrc4) {
        iryrc4 = sryrc4;
    }

    private double iryro = 0; 
    
    //RyR  probability of state O
    public synchronized void setiryro(double sryro) {
        iryro = sryro;
    }

    private double iryri1 = 0.082; 
    
    //RyR  probability of state I1
    public synchronized void setiryri1(double sryri1) {
        iryri1 = sryri1;
    }

    private double iryri2 = 0.002; 
    
    //RyR  probability of state I2
    public synchronized void setiryri2(double sryri2) {
        iryri2 = sryri2;
    }

    private double iryri3 = 0.001; 
    
    //RyR  probability of state I3
    public synchronized void setiryri3(double sryri3) {
        iryri3 = sryri3;
    }

    private double iryri4 = 0; 
    
    //RyR  probability of state I4
    public synchronized void setiryri4(double sryri4) {
        iryri4 = sryri4;
    }

    private double iryri5 = 0; 
    
    //RyR  probability of state I5
    public synchronized void setiryri5(double sryri5) {
        iryri5 = sryri5;
    }

    private double gradedrel; //ICaL-dependent function for determining graded response of Irel

    public double getgradedrel() {
        return gradedrel;
    }

    private double vgainofrel; //Function describing voltage dependence of  Irel gain

    public double getvgainofrel() {
        return vgainofrel;
    }

    private double grel; //Function describing Irel conductance

    public double getgrel() {
        return grel;
    }

    private double scgirel = 1;

    // Scale factor for RyR Current
    public synchronized void setscgirel(double sscgirel) {
        scgirel = sscgirel;
    }

    public double getscgirel() {
        return scgirel;
    }

    /* Current through T-type Ca Channel */
    private double icat; // Ca current through T-type Ca channel (uA/uF)

    public synchronized void seticat(double sicat) {
        icat = sicat;
    }

    public double geticat() {
        return icat;
    }

    private double gcat; // Max. Conductance of the T-type Ca channel (mS/uF)

    public synchronized void setgcat(double sgcat) {
        gcat = sgcat;
    }

    public double getgcat() {
        return gcat;
    }

    private double eca; // Reversal Potential of the T-type Ca channel (mV)

    public synchronized void seteca(double seca) {
        eca = seca;
    }

    public double geteca() {
        return eca;
    }

    private double ab; // Ca alpha-b rate constant (ms^-1)

    public double getab() {
        return ab;
    }

    private double bb; // Ca beta-b rate constant (ms^-1)

    public double getbb() {
        return bb;
    }

    private double b; // Voltage dependent activation gate b

    public double getb() {
        return b;
    }

    private double bss; // Steady-state value of activation gate b

    public double getbss() {
        return bss;
    }

    private double taub; // Time constant of gate b (ms^-1)

    public double gettaub() {
        return taub;
    }

    private double ag; // Ca alpha-g rate constant (ms^-1)

    public double getag() {
        return ag;
    }

    private double bg; // Ca beta-g rate constant (ms^-1)

    public double getbg() {
        return bg;
    }

    private double g; // Voltage dependent inactivation gate g

    public double getg() {
        return g;
    }

    private double gss; // Steady-state value of inactivation gate g

    public double getgss() {
        return gss;
    }

    private double taug; // Time constant of gate g (ms^-1)

    public double gettaug() {
        return taug;
    }

    private double scgcat = 1;

    // Scaling factor for Max. Conductance of the T-type Ca channel
    public synchronized void setscgcat(double sscgcat) {
        scgcat = sscgcat;
    }

    public double getscgcat() {
        return scgcat;
    }

    /* Rapidly Activating Potassium Current */
    private double ikr; // Rapidly Activating K Current (uA/uF)

    public synchronized void setikr(double sikr) {
        ikr = sikr;
    }

    public double getikr() {
        return ikr;
    }

    private double gkr;

    // Channel Conductance of Rapidly Activating K Current (mS/uF)
    public synchronized void setgkr(double sgkr) {
        gkr = sgkr;
    }

    public double getgkr() {
        return gkr;
    }

    private double ekr;

    // Reversal Potential of Rapidly Activating K Current (mV)
    public synchronized void setekr(double sekr) {
        ekr = sekr;
    }

    public double getekr() {
        return ekr;
    }

    private double axr; // K alpha-xr rate constant (ms^-1)

    public double getaxr() {
        return axr;
    }

    private double bxr; // K beta-xr rate constant (ms^-1)

    public double getbxr() {
        return bxr;
    }

    private double xr; // Rapidly Activating K time-dependent activation

    public double getxr() {
        return xr;
    }

    private double xrss; // Steady-state value of inactivation gate xr

    public double getxrss() {
        return xrss;
    }

    private double tauxr; // Time constant of gate xr (ms^-1)

    public double gettauxr() {
        return tauxr;
    }

    private double r; // K time-independent inactivation r

    public double getr() {
        return r;
    }

    private double scgkr = 1;

    // Scale factor for Channel Conductance of Rapidly Activating K Current
    public synchronized void setscgkr(double sscgkr) {
        scgkr = sscgkr;
    }

    public double getscgkr() {
        return scgkr;
    }

    /* Slowly Activating Potassium Current */
    private double iks; // Slowly Activating K Current (uA/uF)

    public synchronized void setiks(double siks) {
        iks = siks;
    }

    public double getiks() {
        return iks;
    }

    private double gks;

    // Channel Conductance of Slowly Activating K Current (mS/uF)
    public synchronized void setgks(double sgks) {
        gks = sgks;
    }

    public double getgks() {
        return gks;
    }

    private double eks;

    // Reversal Potential of Slowly Activating K Current (mV)
    public synchronized void seteks(double seks) {
        eks = seks;
    }

    public double geteks() {
        return eks;
    }

    private double axs1; // K alpha-xs1 rate constant (ms^-1)

    public double getaxs1() {
        return axs1;
    }

    private double axs2; // K alpha-xs2 rate constant (ms^-1)

    public double getaxs2() {
        return axs2;
    }

    private double bxs1; // K beta-xs1 rate constant (ms^-1)

    public double getbxs1() {
        return bxs1;
    }

    private double bxs2; // K beta-xs2 rate constant (ms^-1)

    public double getbxs2() {
        return bxs2;
    }

    private double xs1; // Slowly Activating K time-dependent activation 1

    public double getxs1() {
        return xs1;
    }

    private double xs2; // Slowly Activating K time-dependent activation 2

    public double getxs2() {
        return xs2;
    }

    private double xsss; // Steady-state value of inactivation gate xs

    public double getxsss() {
        return xsss;
    }

    private double tauxs1; // Time constant 1 of gate xs (ms^-1)

    public double gettauxs1() {
        return tauxs1;
    }

    private double tauxs2; // Time constant 2 of gate xs (ms^-1)

    public double gettauxs2() {
        return tauxs2;
    }

    private double prnak = 0.01833; // Na/K Permeability Ratio

    public synchronized void setprnak(double sprnak) {
        prnak = sprnak;
    }

    public double getprnak() {
        return prnak;
    }

    private double scgks = 1;

    // Scaling factor for Channel Conductance of Slowly Activating K Current
    public synchronized void setscgks(double sscgks) {
        scgks = sscgks;
    }

    public double getscgks() {
        return scgks;
    }

    /* Potassium Current (time-independent) */
    private double iki; // Time-independent K current (uA/uF)

    public synchronized void setiki(double siki) {
        iki = siki;
    }

    public double getiki() {
        return iki;
    }

    private double gki;

    // Channel Conductance of Time Independent K Current (mS/uF)
    public synchronized void setgki(double sgki) {
        gki = sgki;
    }

    public double getgki() {
        return gki;
    }

    private double eki;

    // Reversal Potential of Time Independent K Current (mV)
    public synchronized void seteki(double seki) {
        eki = seki;
    }

    public double geteki() {
        return eki;
    }

    private double aki; // K alpha-ki rate constant (ms^-1)

    public double getaki() {
        return aki;
    }

    private double bki; // K beta-ki rate constant (ms^-1)

    public double getbki() {
        return bki;
    }

    private double kin; // K inactivation

    public double getkin() {
        return kin;
    }

    private double scgki = 1;

    // Scaling factor for Channel Conductance of Time Independent K Current
    public synchronized void setscgki(double sscgki) {
        scgki = sscgki;
    }

    public double getscgki() {
        return scgki;
    }

    /* Plateau Potassium Current */
    private double ikp; // Plateau K current (uA/uF)

    public synchronized void setikp(double sikp) {
        ikp = sikp;
    }

    public double getikp() {
        return ikp;
    }

    private double gkp; // Channel Conductance of Plateau K Current (mS/uF)

    public synchronized void setgkp(double sgkp) {
        gkp = sgkp;
    }

    public double getgkp() {
        return gkp;
    }

    private double ekp; // Reversal Potential of Plateau K Current (mV)

    public synchronized void setekp(double sekp) {
        ekp = sekp;
    }

    public double getekp() {
        return ekp;
    }

    private double kp; // K plateau factor

    public double getkp() {
        return kp;
    }

    private double scgkp = 1;

    // Scaling factor for Channel Conductance of Plateau K Current
    public synchronized void setscgkp(double sscgkp) {
        scgkp = sscgkp;
    }

    public double getscgkp() {
        return scgkp;
    }

    /* Sodium-Calcium Exchanger */
    private double inaca; // NaCa exchanger current (uA/uF)

    public synchronized void setinaca(double sinaca) {
        inaca = sinaca;
    }

    public double getinaca() {
        return inaca;
    }

    private double inacass; // NaCa exchanger current in subspace (uA/uF)

    public synchronized void setinacass(double sinacass) {
        inacass = sinacass;
    }

    public double getinacass() {
        return inacass;
    }

    private double knaca = 1; // Scaling factor for inaca (uA/uF) (2000 in LRdII)

    public synchronized void setknaca(double sknaca) {
        knaca = sknaca;
    }

    public double getknaca() {
        return knaca;
    }

    private double kmnae = 87.5;

    // Half-saturation concentration of NaCa exhanger Na (mM)
    public synchronized void setkmnae(double skmnae) {
        kmnae = skmnae;
    }

    public double getkmnae() {
        return kmnae;
    }

    private double kmcae = 1.38;

    // Half-saturation concentration of NaCa exhanger Ca (mM)
    public synchronized void setkmcae(double skmcae) {
        kmcae = skmcae;
    }

    public double getkmcae() {
        return kmcae;
    }

    private double ksat = 0.1;

    // Saturation factor of INaCa at very negative potentials
    public synchronized void setksat(double sksat) {
        ksat = sksat;
    }

    public double getksat() {
        return ksat;
    }

    private double nu = 0.15;//0.35 in LRdII

    // Position of energy barrier controlling voltage dependance of inaca
    public synchronized void setnu(double snu) {
        nu = snu;
    }

    public double getnu() {
        return nu;
    }

    /* Sodium-Potassium Pump */
    private double inak; // NaK pump current (uA/uF)

    public synchronized void setinak(double sinak) {
        inak = sinak;
    }

    public double getinak() {
        return inak;
    }

    private double fnak; // Voltage-dependence parameter of inak

    public double getfnak() {
        return fnak;
    }

    private double sigma; // [Na]o dependence factor of fnak

    public double getsigma() {
        return sigma;
    }

    private double ibarnak = 2.25; // Max. current through Na-K pump (uA/uF) 1.5 in LRdII

    public synchronized void setibarnak(double sibarnak) {
        ibarnak = sibarnak;
    }

    public double getibarnak() {
        return ibarnak;
    }

    private double kmnai = 10;

    // Half-saturation concentration of NaK pump Na (mM)
    public synchronized void setkmnai(double skmnai) {
        kmnai = skmnai;
    }

    public double getkmnai() {
        return kmnai;
    }

    private double kmko = 1.5;

    // Half-saturation concentration of NaK pump K (mM)
    public synchronized void setkmko(double skmko) {
        kmko = skmko;
    }

    public double getkmko() {
        return kmko;
    }

    /* Nonspecific Ca-activated Current */
    private double insna; // Non-specific Na current (uA/uF)

    public synchronized void setinsna(double sinsna) {
        insna = sinsna;
    }

    public double getinsna() {
        return insna;
    }

    private double insk; // Non-specific K current (uA/uF)

    public synchronized void setinsk(double sinsk) {
        insk = sinsk;
    }

    public double getinsk() {
        return insk;
    }

    private double ibarnsna; // Max. Na current through NSCa channel (uA/uF)

    public double getibarnsna() {
        return ibarnsna;
    }

    private double ibarnsk; // Max. K current through NSCa channel (uA/uF)

    public double getibarnsk() {
        return ibarnsk;
    }

    private double pnsca = 0.000000175;

    // Permeability of membrane to Ca (cm/s)
    public synchronized void setpnsca(double spnsca) {
        pnsca = spnsca;
    }

    public double getpnsca() {
        return pnsca;
    }

    private double kmnsca = 0.0012;

    // Half-saturation concentration of NSCa channel (mM)
    public synchronized void setkmnsca(double skmnsca) {
        kmnsca = skmnsca;
    }

    public double getkmnsca() {
        return kmnsca;
    }

    /* ATP-sensitive Potassium Current */
    protected double mgi = 1.0; // Intracellular Mg Concentration (mM)

    public synchronized void setmgi(double smgi) {
        mgi = smgi;
    }

    public double getmgi() {
        return mgi;
    }

    protected double atpi = 6800; // Intracellular ATP Concentration (uM)

    public synchronized void setatpi(double satpi) {
        atpi = satpi;
    }

    public double getatpi() {
        return atpi;
    }

    protected double adpi = 15; // Intracellular ADP Concentration (uM)

    public synchronized void setadpi(double sadpi) {
        adpi = sadpi;
    }

    public double getadpi() {
        return adpi;
    }

    private double ikatp; // ATP-sensitive Potassium Current (uA/uF)

    public synchronized void setikatp(double sikatp) {
        ikatp = sikatp;
    }

    public double getikatp() {
        return ikatp;
    }

    private double fatp; // Fraction of open KATP channels

    public double getfatp() {
        return fatp;
    }

    private double gbarkatp; // IKATP Conductance

    public double getgbarkatp() {
        return gbarkatp;
    }

    private double scgbarkatp = 0; // Scaling factor for IKATP Conductance

    public synchronized void setscgbarkatp(double sscgbarkatp) {
        scgbarkatp = sscgbarkatp;
    }

    public double getscgbarkatp() {
        return scgbarkatp;
    }

    /* Sarcolemmal Ca Pump */
    private double ipca; // Sarcolemmal Ca pump current (uA/uF)

    public synchronized void setipca(double sipca) {
        ipca = sipca;
    }

    public double getipca() {
        return ipca;
    }

    private double ibarpca = 1;

    // Max. Ca current through sarcolemmal Ca pump (uA/uF), 1.15 in LRdII
    public synchronized void setibarpca(double sibarpca) {
        ibarpca = sibarpca;
    }

    public double getibarpca() {
        return ibarpca;
    }

    private double kmpca = 0.0005;

    // Half-saturation concentration of sarcolemmal Ca pump (mM)
    public synchronized void setkmpca(double skmpca) {
        kmpca = skmpca;
    }

    public double getkmpca() {
        return kmpca;
    }

    /* Ca Background Current */
    private double icab; // Ca background current (uA/uF)

    public synchronized void seticab(double sicab) {
        icab = sicab;
    }

    public double geticab() {
        return icab;
    }

    private double gcab; // Max. conductance of Ca background (mS/uF)

    public double getgcab() {
        return gcab;
    }

    private double ecan; // Nernst potential for Ca (mV)

    public synchronized void setecan(double secan) {
        ecan = secan;
    }

    public double getecan() {
        return ecan;
    }

    private double scgcab = 1;

    // Scaling factor for Max. conductance of Ca background
    public synchronized void setscgcab(double sscgcab) {
        scgcab = sscgcab;
    }

    public double getscgcab() {
        return scgcab;
    }

    /* Na Background Current */
    private double inab; // Na background current (uA/uF)

    public synchronized void setinab(double sinab) {
        inab = sinab;
    }

    public double getinab() {
        return inab;
    }

    private double gnab; // Max. conductance of Na background (mS/uF)

    public double getgnab() {
        return gnab;
    }

    private double enan; // Nernst potential for Na (mV)

    public synchronized void setenan(double senan) {
        enan = senan;
    }

    public double getenan() {
        return enan;
    }

    private double scgnab = 1;

    // Scaling factor for Max. conductance of Na background
    public synchronized void setscgnab(double sscgnab) {
        scgnab = sscgnab;
    }

    public double getscgnab() {
        return scgnab;
    }

    /* ** CONSTRUCTOR ** */
    public FaberRudyModel_I_2007() {
        super();
        getCellGeometry();
    }

    public void setCurrents() {
        currentINa();
        currentICal();
        currentIrel();
        currentICat();
        currentIKr();
        currentIKs();
        currentIKi();
        currentIKp();
        currentINaCa();
        currentINaK();
        /* This current is only used when modeling pathological conditions */
        //currentINasCa ();
        /* This current is only used when modeling pathological conditions */
        //currentIKatp();
        currentIpCa();
        currentICab();
        currentINab();
        currentIt();
    }

    public void setConcentrations() {
        SRfluxes();
        CabuffMyo();
        CabuffJsr();
        concentrationMyo();
        concentrationNsr();
        concentrationJsr();
        concentrationSs();
    }

    public void doSimulate() {
        setCurrents();
        if (t >= 0.005) setConcentrations();

        vnew = v - it * dt;
        dvdtnew = (vnew - v) / dt;

        //getCaRelease();

        t = t + dt;
        v = vnew;
        dvdt = dvdtnew;
    }

    /** ***************************************************** */
    /**
     * Cell Geometry
     */
    protected void getCellGeometry() {
        vcell = 1000 * Math.PI * a * a * l;
        ageo = 2 * Math.PI * a * a + 2 * Math.PI * a * l;
        acap = ageo * 2;
        vmyo = vcell * 0.68;
        vmito = vcell * 0.24;
        vsr = vcell * 0.06;
        vnsr = vcell * 0.0552;
        vjsr = vcell * 0.0048;
        //vcleft = vcell * 0.12 / 0.88;
        vss = vcell * 0.02;
    }

    /**
     * This is where CICR is calculated, following dv/dtmax cai is measured, 2
     * ms later cai is measured again, the difference of the two determines the
     * amount of SR release
     */
    //protected void getCaRelease() {
    //    if ((v > -35.0) && (dvdtnew < dvdt) && (dvdtmax == false)) {
    //        caitzero = cai + cmdn + trpn + dcai + csqn + jsr + nsr;
    //        dvdtmax = true;
    //        tcount = 0;
    //    }

    //    if ((tcount >= 2) && (tcount < (2 + dt))) {
    //        caiatwo = cai + cmdn + trpn + dcai + csqn + jsr + nsr;
    //        dcaitwo = caiatwo - caitzero;
    //        if (dcaitwo >= dcaith) {
    //            grelbarcicr = gmaxrel * (dcaitwo - 0.00018)
    //                    / (kmrel + dcaitwo - 0.00018);
    //            tcicr = 0;
    //        } else if (dcaitwo < 0) {
    //            grelbarcicr = 0;
    //        } else if ((dcaitwo >= 0) && (dcaitwo < dcaith)) {
    //            dcai = dcaitwo * 1000.0;
    //            grelbarcicr = gmaxrel
    //                    * (7.5 * (Math.pow(dcai, 3.0)) - (Math.pow(dcai, 2.0)) + 0.1 * (dcai));
    //        }
    //    }

    //    tcount = tcount + dt;

        /* This is the function that determines if a spontaneous release occurs */

    //    if ((csqn >= csqnth) && (tjsrol > 50.0)) {
    //        grelbarjsrol = 4.0;
    //        tjsrol = 0.0;
    //    }
    //}

    /* Functions that describe the currents begin here */

    /**
     * Calculates Fast Na Current
     */
    protected void currentINa() {
        gna = 16.0 * scgna;
        ena = ((R * temp) / frdy) * Math.log(nao / nai);

        if (v != -47.13) {
            am = 0.32 * (v + 47.13) / (1 - Math.exp(-0.1 * (v + 47.13)));
        } else {
            am = 3.2;
        }
        bm = 0.08 * Math.exp(-v / 11.0);

        if (v < -40.0) {
            ah = 0.135 * Math.exp((80.0 + v) / -6.8);
            bh = 3.56 * Math.exp(0.079 * v) + 310000.0 * Math.exp(0.35 * v);
            aj = (-127140.0 * Math.exp(0.2444 * v) - 0.00003474 * Math
                    .exp(-0.04391 * v))
                    * (v + 37.78) / (1.0 + Math.exp(0.311 * (v + 79.23)));
            bj = (0.1212 * Math.exp(-0.01052 * v))
                    / (1.0 + Math.exp(-0.1378 * (v + 40.14)));
        } else {
            ah = 0.0;
            bh = 1.0 / (0.13 * (1.0 + Math.exp((v + 10.66) / -11.1)));
            aj = 0.0;
            bj = (0.3 * Math.exp(-0.0000002535 * v))
                    / (1.0 + Math.exp(-0.1 * (v + 32.0)));
        }

        if (t == 0) // Steady State Values
        {
            h = ah / (ah + bh);
            j = aj / (aj + bj);
            m = am / (am + bm);
        } else {
            h = ah / (ah + bh) - ((ah / (ah + bh)) - h)
                    * Math.exp(-dt / (1 / (ah + bh)));
            j = aj / (aj + bj) - ((aj / (aj + bj)) - j)
                    * Math.exp(-dt / (1 / (aj + bj)));
            m = am / (am + bm) - ((am / (am + bm)) - m)
                    * Math.exp(-dt / (1 / (am + bm)));
        }

        ina = gna * m * m * m * h * j * (v - ena);
    }

    /**
     * Calculates Currents through L-Type Ca Channel
     */
    protected void currentICal() {
        //dss = 1.0 / (1.0 + Math.exp(-(v + 10.0) / 6.24));
        //if (v != -10.0) {
        //    taud = dss * (1.0 - Math.exp(-(v + 10.0) / 6.24))
        //            / (0.035 * (v + 10.0));
        //} else {
        //    taud = dss * (1 / 6.24) / 0.035;
        //}
        //ad = dss / taud;
        //bd = (1 - dss) / taud;

        //fss = (1.0 / (1.0 + Math.exp((v + 32.0) / 8.0)))
        //        + (0.6 / (1.0 + Math.exp((50.0 - v) / 20.0)));
        //tauf = 1.0 / (0.0197 * Math.exp(-Math.pow(0.0337 * (v + 10.0), 2.0)) + 0.02);
        //af = fss / tauf;
        //bf = (1 - fss) / tauf;
        icala = 0.925 * Math.exp(v/30);
        icalb = 0.39 * Math.exp(-v/40);
        icalgf = 0.245 * Math.exp(v/10);
        icalgs = 0.005 * Math.exp(-v/40);
        icalff = 0.02 * Math.exp(v/500);
        icalfs = 0.03 * Math.exp(-v/280);
        icallf = 0.035 * Math.exp(-v/300);
        icalls = 0.0011 * Math.exp(v/500);
        icalof = (4 * icalb * icallf * icalgf) / (icala * icalff);
        icalos = (4 * icalb * icalls * icalgs) / (icala * icalfs);
        icalosf = (icalls * icalff) / icallf;
        icalofs = icalfs;
        icald = 4 / (1 + 1/cass);
        icalt = 0.01; 

        if (t == 0) // Steady State Values
        {
            int iindex = 0;
            int ilength = 30000;
            for (iindex = 0; iindex < ilength; iindex++){
                   double newicalc0 = icalc0 + (icalc1 * icalb + icalc0ca * icalt - icalc0 * (icala + icald)) * dt;
            double newicalc1 = icalc1 + (icalc0 * icala + icalc2 * 2 * icalb + icalc1ca * icalt - icalc1 
                    * (3 * icala + icalb + icald)) * dt;
            double newicalc2 = icalc2 + (icalc1 * 3 * icala + icalc3 * 3 * icalb + icalc2ca * icalt - icalc2 
                    * (2 * icala + 2 * icalb + icald)) * dt;
            double newicalc3 = icalc3 + (icalc2 * 2 * icala + icalo * 4 * icalb + icalif * icalof + icalis * icalos
                    + icalc3ca * icalt - icalc3 * (icala + 3 * icalb + icalgf + icalgs + icald)) * dt;
            double newicalo = icalo + (icalc3 * icala + icalif * icallf + icalis * icalls + icalica * icalt
                    - icalo * (4 * icalb + icalff + icalfs + icald)) * dt;
            double newicalif = icalif + (icalc3 * icalgf + icalo * icalff + icalis * icalosf + icalifca * icalt
                    - icalif * (icalof + icallf + icalofs + icald)) * dt;
            double newicalis = icalis + (icalc3 * icalgs + icalo * icalfs + icalif * icalofs + icalisca * icalt
                    - icalis * (icalos + icalls + icalosf + icald)) * dt;
            double newicalc0ca = icalc0ca + (icalc1ca * icalb + icalc0 * icald - icalc0ca * (4 * icala + icalt)) * dt;
            double newicalc1ca = icalc1ca + (icalc0ca * 4 * icala + icalc2ca * 2 * icalb + icalc1 * icald
                    - icalc1ca * (3 * icala + icalb + icalt)) * dt;
            double newicalc2ca = icalc2ca + (icalc1ca * 3 * icala + icalc3ca * 3 * icalb + icalc2 * icald
                    - icalc2ca * (2 * icala + 2 * icalb + icalt)) * dt;
            double newicalc3ca = icalc3ca + (icalc2ca * 2 * icala + icalica * 4 * icalb + icalifca * icalof
                    + icalisca * icalos +icalc3 * icald - icalc3ca * (icala + 3 * icalb + icalgf + icalgs + icalt)) * dt;
            double newicalica = icalica + (icalc3ca * icala + icalifca * icallf + icalisca * icalls + icalo * icald
                    - icalica * (4 * icalb + icalff + icalfs + icalt)) * dt;
            double newicalifca = icalifca + (icalc3ca * icalgf + icalica * icalff + icalisca * icalosf + icalif * icald
                    - icalifca * (icalof + icallf + icalofs + icalt)) * dt;
            double newicalisca = icalisca + (icalc3ca * icalgs + icalica * icalfs + icalifca * icalofs + icalis * icald
                    - icalisca * (icalos + icalls + icalosf + icalt)) * dt;
            icalc0 = newicalc0;
            icalc1 = newicalc1;
            icalc2 = newicalc2;
            icalc3 = newicalc3;
            icalo = newicalo;
            icalif = newicalif;
            icalis = newicalis;
            icalc0ca = newicalc0ca;
            icalc1ca = newicalc1ca;
            icalc2ca = newicalc2ca;
            icalc3ca = newicalc3ca;
            icalica = newicalica;
            icalifca = newicalifca;
            icalisca = newicalisca;
            }
        //    d = ad / (ad + bd);
        //    f = af / (af + bf);
        } else {
            double newicalc0 = icalc0 + (icalc1 * icalb + icalc0ca * icalt - icalc0 * (icala + icald)) * dt;
            double newicalc1 = icalc1 + (icalc0 * icala + icalc2 * 2 * icalb + icalc1ca * icalt - icalc1 
                    * (3 * icala + icalb + icald)) * dt;
            double newicalc2 = icalc2 + (icalc1 * 3 * icala + icalc3 * 3 * icalb + icalc2ca * icalt - icalc2 
                    * (2 * icala + 2 * icalb + icald)) * dt;
            double newicalc3 = icalc3 + (icalc2 * 2 * icala + icalo * 4 * icalb + icalif * icalof + icalis * icalos
                    + icalc3ca * icalt - icalc3 * (icala + 3 * icalb + icalgf + icalgs + icald)) * dt;
            double newicalo = icalo + (icalc3 * icala + icalif * icallf + icalis * icalls + icalica * icalt
                    - icalo * (4 * icalb + icalff + icalfs + icald)) * dt;
            double newicalif = icalif + (icalc3 * icalgf + icalo * icalff + icalis * icalosf + icalifca * icalt
                    - icalif * (icalof + icallf + icalofs + icald)) * dt;
            double newicalis = icalis + (icalc3 * icalgs + icalo * icalfs + icalif * icalofs + icalisca * icalt
                    - icalis * (icalos + icalls + icalosf + icald)) * dt;
            double newicalc0ca = icalc0ca + (icalc1ca * icalb + icalc0 * icald - icalc0ca * (4 * icala + icalt)) * dt;
            double newicalc1ca = icalc1ca + (icalc0ca * 4 * icala + icalc2ca * 2 * icalb + icalc1 * icald
                    - icalc1ca * (3 * icala + icalb + icalt)) * dt;
            double newicalc2ca = icalc2ca + (icalc1ca * 3 * icala + icalc3ca * 3 * icalb + icalc2 * icald
                    - icalc2ca * (2 * icala + 2 * icalb + icalt)) * dt;
            double newicalc3ca = icalc3ca + (icalc2ca * 2 * icala + icalica * 4 * icalb + icalifca * icalof
                    + icalisca * icalos +icalc3 * icald - icalc3ca * (icala + 3 * icalb + icalgf + icalgs + icalt)) * dt;
            double newicalica = icalica + (icalc3ca * icala + icalifca * icallf + icalisca * icalls + icalo * icald
                    - icalica * (4 * icalb + icalff + icalfs + icalt)) * dt;
            double newicalifca = icalifca + (icalc3ca * icalgf + icalica * icalff + icalisca * icalosf + icalif * icald
                    - icalifca * (icalof + icallf + icalofs + icalt)) * dt;
            double newicalisca = icalisca + (icalc3ca * icalgs + icalica * icalfs + icalifca * icalofs + icalis * icald
                    - icalisca * (icalos + icalls + icalosf + icalt)) * dt;
            icalc0 = newicalc0;
            icalc1 = newicalc1;
            icalc2 = newicalc2;
            icalc3 = newicalc3;
            icalo = newicalo;
            icalif = newicalif;
            icalis = newicalis;
            icalc0ca = newicalc0ca;
            icalc1ca = newicalc1ca;
            icalc2ca = newicalc2ca;
            icalc3ca = newicalc3ca;
            icalica = newicalica;
            icalifca = newicalifca;
            icalisca = newicalisca;

        //    d = ad / (ad + bd) - ((ad / (ad + bd)) - d)
        //            * Math.exp(-dt / (1 / (ad + bd)));
        //    f = af / (af + bf) - ((af / (af + bf)) - f)
        //            * Math.exp(-dt / (1 / (af + bf)));
        }

        if (v != 0) {
            ibarca = pca
                    * zca
                    * zca
                    * ((v * frdy * frdy) / (R * temp))
                    * ((gacai * cass * Math.exp((zca * v * frdy) / (R * temp)) - gacao
                            * cao) / (Math.exp((zca * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibarca = pca * 2 * frdy * (gacai * cass - gacao * cao);
        }
        if (v != 0) {
            ibarna = pna
                    * zna
                    * zna
                    * ((v * frdy * frdy) / (R * temp))
                    * ((ganai * nass * Math.exp((zna * v * frdy) / (R * temp)) - ganao
                            * nao) / (Math.exp((zna * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibarna = pna * frdy * (ganai * nass - ganao * nao);
        }
        if (v != 0) {
            ibark = pk
                    * zk
                    * zk
                    * ((v * frdy * frdy) / (R * temp))
                    * ((gaki * kss * Math.exp((zk * v * frdy) / (R * temp)) - gako
                            * ko) / (Math.exp((zk * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibark = pk * frdy * (gaki * kss - gako * ko);
        }

   //     fca = 1.0 / (1.0 + cai / kmca);

        ilca = icalo * ibarca * scgcal;
        ilcana = icalo * ibarna * scgcal;
        ilcak = icalo * ibark * scgcal;

        ilcatot = ilca + ilcana + ilcak;
    }

    /**
     * Calculates Currents through Ryanodyne Receptor Channel
     */
    protected void currentIrel() {
        iryra1 = 1750 * cass;
        iryra2 = 5600 * cass;
        iryrb1 = 5;
        iryrb2 = 2.625;
        iryrb3 = 1;
        iryrb4 = 6.25;
        iryrg1 = 0.4 * cass;
        iryrg2 = 1.2 * cass;
        iryrg3 = 2.8 * cass;
        iryrg4 = 5.2 * cass;
        iryrg5 = 8.4 * cass;
        iryrd1 = 0.01 / (1 + Math.pow((csqnbar * 0.75) / csqn, 9));
        iryrd2 = iryrd1 / 10;
        iryrd3 = iryrd1 / 100;
        iryrd4 = iryrd1 / 1000;
        iryrd5 = iryrd1 / 10000;

        if (t == 0) // Steady State Values
        {
            int kindex = 0;
            int klength = 30000;
            for (kindex = 0; kindex < klength; kindex++){
                   double newiryrc1 = iryrc1 + (iryrc2 * iryrb1 + iryri1 * iryrd1 - iryrc1 * (iryra1 + iryrg1)) * dt;
            double newiryrc2 = iryrc2 + (iryrc1 * iryra1 + iryrc3 * iryrb2 + iryri2 * iryrd2 
                   - iryrc2 * (iryrb1 + iryra2 + iryrg2)) * dt;
            double newiryrc3 = iryrc3 + (iryrc2 * iryra2 + iryrc4 * iryrb3 + iryri3 * iryrd3 
                   - iryrc3 * (iryrb2 + iryra2 + iryrg3)) * dt;
            double newiryrc4 = iryrc4 + (iryrc3 * iryra2 + iryro * iryrb4 + iryri4 * iryrd4 
                   - iryrc4 * (iryrb3 + iryra2 + iryrg4)) * dt;
            double newiryro = iryro + (iryrc4 * iryra2 + iryri5* iryrd5 - iryro * (iryrb4 + iryrg5)) * dt;
            double newiryri1 = iryri1 + (iryri2 * iryrb1 + iryrc1* iryrg1 - iryri1 * (iryra1 + iryrd1)) * dt;
            double newiryri2 = iryri2 + (iryri1 * iryra1 + iryri3 * iryrb2 + iryrc2 * iryrg2 
                   - iryri2 * (iryrb1 + iryra2 + iryrd2)) * dt;
            double newiryri3 = iryri3 + (iryri2 * iryra2 + iryri4 * iryrb3 + iryrc3 * iryrg3 
                   - iryri3 * (iryrb2 + iryra2 + iryrd3)) * dt;
            double newiryri4 = iryri4 + (iryri3 * iryra2 + iryri5 * iryrb4 + iryrc4 * iryrg4 
                   - iryri4 * (iryrb3 + iryra2 + iryrd4)) * dt;
            double newiryri5 = iryri5 + (iryri4 * iryra2 + iryro* iryrg5 - iryri5 * (iryrb4 + iryrd5)) * dt;
            iryrc1 = newiryrc1;
            iryrc2 = newiryrc2;
            iryrc3 = newiryrc3;
            iryrc4 = newiryrc4;
            iryro = newiryro;
            iryri1 = newiryri1;
            iryri2 = newiryri2;
            iryri3 = newiryri3;
            iryri4 = newiryri4;
            iryri5 = newiryri5;
            }
        } else {
            double newiryrc1 = iryrc1 + (iryrc2 * iryrb1 + iryri1 * iryrd1 - iryrc1 * (iryra1 + iryrg1)) * dt;
            double newiryrc2 = iryrc2 + (iryrc1 * iryra1 + iryrc3 * iryrb2 + iryri2 * iryrd2 
                   - iryrc2 * (iryrb1 + iryra2 + iryrg2)) * dt;
            double newiryrc3 = iryrc3 + (iryrc2 * iryra2 + iryrc4 * iryrb3 + iryri3 * iryrd3 
                   - iryrc3 * (iryrb2 + iryra2 + iryrg3)) * dt;
            double newiryrc4 = iryrc4 + (iryrc3 * iryra2 + iryro * iryrb4 + iryri4 * iryrd4 
                   - iryrc4 * (iryrb3 + iryra2 + iryrg4)) * dt;
            double newiryro = iryro + (iryrc4 * iryra2 + iryri5* iryrd5 - iryro * (iryrb4 + iryrg5)) * dt;
            double newiryri1 = iryri1 + (iryri2 * iryrb1 + iryrc1* iryrg1 - iryri1 * (iryra1 + iryrd1)) * dt;
            double newiryri2 = iryri2 + (iryri1 * iryra1 + iryri3 * iryrb2 + iryrc2 * iryrg2 
                   - iryri2 * (iryrb1 + iryra2 + iryrd2)) * dt;
            double newiryri3 = iryri3 + (iryri2 * iryra2 + iryri4 * iryrb3 + iryrc3 * iryrg3 
                   - iryri3 * (iryrb2 + iryra2 + iryrd3)) * dt;
            double newiryri4 = iryri4 + (iryri3 * iryra2 + iryri5 * iryrb4 + iryrc4 * iryrg4 
                   - iryri4 * (iryrb3 + iryra2 + iryrd4)) * dt;
            double newiryri5 = iryri5 + (iryri4 * iryra2 + iryro* iryrg5 - iryri5 * (iryrb4 + iryrd5)) * dt;
            iryrc1 = newiryrc1;
            iryrc2 = newiryrc2;
            iryrc3 = newiryrc3;
            iryrc4 = newiryrc4;
            iryro = newiryro;
            iryri1 = newiryri1;
            iryri2 = newiryri2;
            iryri3 = newiryri3;
            iryri4 = newiryri4;
            iryri5 = newiryri5;
            
            gradedrel = 1 / (1 + Math.exp((20 + ilca) / 6)) - 0.034445; 
            vgainofrel = 1 + Math.exp((0.015 * ibarca + 1.25) / 0.75);
            grel = (250 * gradedrel) / vgainofrel;

            irel = scgirel * grel * iryro * (jsr - cass);
            }
    }
    /**
     * Calculates Currents through T-Type Ca Channel
     */
    protected void currentICat() {
        bss = 1.0 / (1.0 + Math.exp(-(v + 14.0) / 10.8));
        taub = 3.7 + 6.1 / (1.0 + Math.exp((v + 25.0) / 4.5));
        ab = bss / taub;
        bb = (1 - bss) / taub;

        gss = 1.0 / (1.0 + Math.exp((v + 60.0) / 5.6));
        if (v <= 0)
            taug = -0.875 * v + 12.0;
        else
            taug = 12.0;
        ag = gss / taug;
        bg = (1 - gss) / taug;

        if (t == 0) // Steady State Values
        {
            b = ab / (ab + bb);
            g = ag / (ag + bg);
        } else {
            b = ab / (ab + bb) - ((ab / (ab + bb)) - b)
                    * Math.exp(-dt / (1 / (ab + bb)));
            g = ag / (ag + bg) - ((ag / (ag + bg)) - g)
                    * Math.exp(-dt / (1 / (ag + bg)));
        }

        gcat = 0.005 * scgcat;
        eca = ((R * temp) / (2.0 * frdy)) * Math.log(cao / cai);

        icat = gcat * b * b * g * (v - eca);
    }

    /**
     * Calculates Rapidly Activating K Current
     */
    protected void currentIKr() {
        gkr = 0.031368 * Math.sqrt(ko / 5.4);
	if (scgkr != 1) {
	  gkr *= scgkr; // gkr = gkr * scgkr;
	}
        ekr = ((R * temp) / frdy) * Math.log(ko / ki);

        xrss = 1.0 / (1.0 + Math.exp(-(v + 10.085) / 4.25));
        if (v != -14.2 || v != -38.9) {
            tauxr = 1.0 / (0.00138 * (v + 14.2)
                    / (1.0 - Math.exp(-0.123 * (v + 14.2))) + 0.00061
                    * (v + 38.9) / (Math.exp(0.145 * (v + 38.9)) - 1.0));
        } else {
            if (v == -14.2) {
                tauxr = 1.0 / (0.00138 / 0.123 + 0.00061 * (v + 38.9)
                        / (Math.exp(0.145 * (v + 38.9)) - 1.0));
            } else {
                if (v == -38.9) {
                    tauxr = 1.0 / (0.00138 * (v + 14.2)
                            / (1.0 - Math.exp(-0.123 * (v + 14.2))) + 0.00061 / 0.145);
                }
            }
        }

        axr = xrss / tauxr;
        bxr = (1 - xrss) / tauxr;

        if (t == 0) // Steady State Values
        {
            xr = axr / (axr + bxr);
        } else {
            xr = axr / (axr + bxr) - ((axr / (axr + bxr)) - xr)
                    * Math.exp(-dt / (1 / (axr + bxr)));
        }

        r = 1.0 / (1.0 + Math.exp((v + 9.0) / 18.4));//22.4 in LRdII

        ikr = gkr * xr * r * (v - ekr);
    }

    /**
     * Calculates Slowly Activating K Current
     */
    protected void currentIKs() {
        gks = 0.3031 * (1.0 + 0.6 / (1.0 + Math.pow(0.000038 / cai, 1.4)));
	if (scgks != 1) {
	  gks *= scgks;
	}
        eks = ((R * temp) / frdy)
                * Math.log((ko + prnak * nao) / (ki + prnak * nai));

        xsss = 1.0 / (1.0 + Math.exp(-(v - 1.5) / 16.7));

        if (v != -30.0) {
            tauxs1 = 1.0 / (0.0000719 * (v + 30.0)
                    / (1.0 - Math.exp(-0.148 * (v + 30.0))) + 0.000131
                    * (v + 30.0) / (Math.exp(0.0687 * (v + 30.0)) - 1.0));
        } else {
            tauxs1 = 1.0 / 0.0000719 / 0.148;
        }
        tauxs2 = tauxs1 * 4;

        axs1 = xsss / tauxs1;
        axs2 = xsss / tauxs2;
        bxs1 = (1 - xsss) / tauxs1;
        bxs2 = (1 - xsss) / tauxs2;

        if (t == 0) // Steady State Values
        {
            xs1 = axs1 / (axs1 + bxs1);
            xs2 = axs2 / (axs2 + bxs2);
        } else {
            xs1 = axs1 / (axs1 + bxs1) - ((axs1 / (axs1 + bxs1)) - xs1)
                    * Math.exp(-dt / (1 / (axs1 + bxs1)));
            xs2 = axs2 / (axs2 + bxs2) - ((axs2 / (axs2 + bxs2)) - xs2)
                    * Math.exp(-dt / (1 / (axs2 + bxs2)));
        }

        iks = gks * xs1 * xs2 * (v - eks);
    }

    /**
     * Calculates Time-Independent K Current IK1
     */
    protected void currentIKi() {
        gki = 0.5625 * Math.sqrt(ko / 5.4);
	if (scgki != 1) {
	  gki *= scgki;
	}
        eki = ((R * temp) / frdy) * Math.log(ko / ki);

        aki = 1.02 / (1.0 + Math.exp(0.2385 * (v - eki - 59.215)));
        bki = (0.49124 * Math.exp(0.08032 * (v - eki + 5.476)) + Math
                .exp(0.06175 * (v - eki - 594.31)))
                / (1.0 + Math.exp(-0.5143 * (v - eki + 4.753)));

        kin = aki / (aki + bki);

        iki = gki * kin * (v - eki);
    }

    /**
     * Calculates Plateau K Current
     */
    protected void currentIKp() {
        gkp = 0.000552 * scgkp;
        ekp = eki;

        kp = 1.0 / (1.0 + Math.exp((7.488 - v) / 5.98));

        ikp = gkp * kp * (v - ekp);
    }

    /**
     * Calculates Na-Ca Exchanger Current
     */
    protected void currentINaCa() {
        //inaca = knaca
        //        * (1.0 / (Math.pow(kmnae, 3.0) + Math.pow(nao, 3.0)))
        //        * (1.0 / (kmcae + cao))
        //        * (1.0 / (1.0 + ksat
        //                * Math.exp(((nu - 1.0) * v * frdy) / (R * temp))))
        //        * (Math.pow(nai, 3.0) * cao
        //                * Math.exp((nu * v * frdy) / (R * temp)) - Math.pow(
        //                nao, 3.0)
        //                * cai * Math.exp(((nu - 1.0) * v * frdy) / (R * temp)));
        inaca = (knaca * 0.8 * 0.000225 * Math.exp(((nu - 1) * v * frdy) / (R * temp))
                * (Math.exp((v * frdy) / (R * temp)) * Math.pow(nai, 3.0) * cao - 1.5 * Math.pow(nao, 3.0) * cai))/
                (1 + 0.0001 * Math.exp(((nu - 1) * v * frdy) / (R * temp))
                * (Math.exp((v * frdy) / (R * temp)) * Math.pow(nai, 3.0) * cao 
                + 1.5 * Math.pow(nao, 3.0) * cai));
        inacass = (knaca * 0.2 * 0.000225 * Math.exp(((nu - 1) * v * frdy) / (R * temp))
                * (Math.exp((v * frdy) / (R * temp)) * Math.pow(nass, 3.0) * cao - 1.5 * Math.pow(nao, 3.0) * cass))/
                (1 + 0.0001 * Math.exp(((nu - 1) * v * frdy) / (R * temp))
                * (Math.exp((v * frdy) / (R * temp)) * Math.pow(nass, 3.0) * cao 
                + 1.5 * Math.pow(nao, 3.0) * cass));
     }

    /**
     * Calculates Na-K Pump Current
     */
    protected void currentINaK() {
        sigma = (Math.exp(nao / 67.3) - 1.0) / 7.0;

        fnak = 1.0 / (1.0 + 0.1245 * Math.exp(- (0.1 * v * frdy) / (R * temp)) + 0.0365
                * sigma * Math.exp(- (v * frdy) / (R * temp)));

        inak = ibarnak * fnak * (1.0 / (1.0 + Math.pow(kmnai / nai, 2.0)))
                * (ko / (ko + kmko));
    }

    /**)
     * Calculates Non-Specific Ca-Activated Current
     */
    protected void currentINasCa() {
        if (v != 0.0) {
            ibarnsna = pnsca
                    * zna
                    * zna
                    * ((v * frdy * frdy) / (R * temp))
                    * ((ganai * nai * Math.exp((zna * v * frdy) / (R * temp)) - ganao
                            * nao) / (Math.exp((zna * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibarnsna = pnsca * frdy * (ganai * nai - ganao * nao);
        }
        if (v != 0.0) {
            ibarnsk = pnsca
                    * zk
                    * zk
                    * ((v * frdy * frdy) / (R * temp))
                    * ((gaki * ki * Math.exp((zk * v * frdy) / (R * temp)) - gako
                            * ko) / (Math.exp((zk * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibarnsk = pnsca * frdy * (gaki * ki - gako * ko);
        }

        insna = ibarnsna / (1.0 + Math.pow(kmnsca / cai, 3.0));
        insk = ibarnsk / (1.0 + Math.pow(kmnsca / cai, 3.0));
    }

    /**
     * Calculates ATP-sensitive Potassium Current
     * <p>
     * Ferrero et al., Circ. Res. 1997;79:208-21.
     */
    protected void currentIKatp() {
        //channel density
        double satp = 3.8;
        //maximum channel probability
        double p0 = 0.91;

        //unitary conductance in the absence of ATP
        //double gamma0 = 35.375 * (Math.pow((ko / 5.4), 0.24));
        //modified to account for the steeper Ko-dependence
        double gamma0 = 35.375 * (Math.pow(ko / 5.4, 0.54));

        //rectification caused by Mg
        double khmg = (0.65 / Math.sqrt(ko + 5.0))
                * Math.exp(-(0.64 * frdy * v) / (R * temp));
        double fm = 1.0 / (1.0 + (mgi / khmg));

        //rectification caused by Na
        double khna = 25.9 * Math.exp(-(0.35 * frdy * v) / (R * temp));
        double fn = 1.0 / (1.0 + Math.pow(nai / khna, 2.0));

        //Temperature effect
        double ft = 1.0; //Math.pow(1.3, ((temp - 36.0) / 10.0));

        //unitary conductance
        double g0 = gamma0 * fm * fn * ft;

        //fraction of open KATP channels
        double kmatp = 35.8 + 17.9 * Math.pow(adpi, 0.256);
        double hatp = 1.3 + 0.74 * Math.exp(-0.09 * adpi);

        fatp = 1.0 / (1.0 + Math.pow(atpi / kmatp, hatp));

        gbarkatp = satp * scgbarkatp;

        //From Shawn, Rudy, Cardiovasc. Res. 1997;35:256-72
        //patp = 1 / (1 + (Math.pow(3 / .25, 2)));
        //ikatp = gbarkatp * patp * (Math.pow((ko / 4), 0.24)) * (v - eki);

        ikatp = gbarkatp * g0 * p0 * fatp * (v - eki);
    }

    /**
     * Calculates Sarcolemmal Ca Pump
     */
    protected void currentIpCa() {
        //ipca = (ibarpca * cai) / (kmpca + cai);
        ipca = (ibarpca * cai) / ((kmpca + cai) * (1 + Math.exp((0.00012 - cai) / 0.00001)));
    }

    /**
     * Calculates Ca Background Current
     */
    protected void currentICab() {
        gcab = 0.003016 * scgcab;
        ecan = ((R * temp) / (2.0 * frdy)) * Math.log(cao / cai);

        icab = gcab * (v - ecan);
    }

    /**
     * Calculates Na Background Current
     */
    protected void currentINab() {
        gnab = 0.0002 * scgnab;
        enan = ena;

        inab = gnab * (v - enan);
    }

    /**
     * Calculates Total Na Current
     */
    protected double getTotalINa() {
        //return getina() + getinab() + getilcana() + getinsna() + 3 * getinak()
        return getina() + getinab() + getilcana() + 3 * getinak()
                + 3 * getinaca() + 3 * getinacass();
    }

    /**
     * Calculates Total K Current
     */
    protected double getTotalIK() {
        return getikr() + getiks() + getiki() + getikp() + getilcak()
                    - 2 * getinak();
               //        + getinsk() + getikatp() - 2 * getinak();
    }

    /**
     * Calculates Total Ca Current
     */
    protected double getTotalICa() {
        return getilca() + geticat() + geticab() + getipca() - 2 * getinaca() - 2 * getinacass();
    }

    /**
     * Calculates Total Current
     *
     * Total sum of currents is calculated here, if the time is between stimtime =
     * 0 and stimtime = tstim, a stimulus is applied
     */
    public void currentIt() {
        naiont = getTotalINa();
        kiont = getTotalIK();
        caiont = getTotalICa();

        if ((t >= tstim) && (t < (tstim + dt))) {
            stimtime = 0;
            tstim = tstim + bcl;
            dvdtmax = false;
        }

        if ((stimtime >= 0) && (stimtime < stl)) {
            it = st + naiont + kiont + caiont;
        } else {
            it = naiont + kiont + caiont;
        }

        stimtime = stimtime + dt;
    }

    /* Functions that calculate intracellular ion concentrations */

    /**
     * Calculates new SR ion fluxes
     */
    protected void SRfluxes() {
        iup = (iupbar * cai) / (cai + kmup);//Ca uptake from myoplasm to SR, mM/ms

 //       kleak = iupbar / nsrbar;
        kleak = 0.00875 / nsrbar; 
        ileak = kleak * nsr;//Ca leak from NSR to myoplasm, mM/ms

        itr = (nsr - jsr) / tautr;//Ca transfer form NSR to JSR, mM/ms    
    }

    /**
     * Calculates new Ca buffer concentrations in myoplasm
     */
    protected void CabuffMyo() {
        trpn = trpnbar * (cai / (cai + kmtrpn));//troponin, Ca buffer in myoplasm
        cmdn = cmdnbar * (cai / (cai + kmcmdn));//calmodulin, Ca buffer in myoplasm  
    }

    /**
     * Calculates new Ca buffer concentrations in JSR
     */
    protected void CabuffJsr() {
       csqn = csqnbar * (jsr / (jsr + kmcsqn));//calsequestrin, Ca buffer in JSR    
    }

    /**
     * Calculates new myoplasmic ion concentrations
     */
    protected void concentrationMyo() {
        bmyo = 1 / (1 + (trpnbar * kmtrpn) / (Math.pow((kmtrpn + cai), 2)) + (cmdnbar * kmcmdn) / 
               (Math.pow((kmcmdn + cai), 2)));
        dcai = -dt * bmyo * ((caiont - ilca + 2 * inacass) * (acap / (vmyo * zca * frdy)) + (iup - ileak) * (vnsr / vmyo)
               - idiffca * (vss / vmyo));
        cai = dcai + cai;
        catotal = trpn + cmdn + cai;
     //    dnai = -dt * (naiont * acap) / (vmyo * zna * frdy);
        dnai = -dt * ((naiont - ilcana - 3 * inacass) * (acap / (vmyo * zna * frdy)) - idiffna * (vss / vmyo));
        nai = dnai + nai;
    //    dki = -dt * (kiont * acap) / (vmyo * zk * frdy);
        dki = -dt * ((kiont - ilcak) * (acap / (vmyo * zk * frdy)) - idiffk * (vss / vmyo));
        ki = dki + ki;
    }

    /**
     * Calculates new NSR Ca ion concentration
     */
    protected void concentrationNsr() {
        dnsr = dt * (iup - ileak - itr * (vjsr / vnsr));
        nsr = nsr + dnsr;
    }

    /**
     * Calculates new JSR Ca ion concentration
     */
    protected void concentrationJsr() {
    //    grelcicr = grelbarcicr * (1.0 - Math.exp(-tcicr / tauon))
    //            * Math.exp(-tcicr / tauoff);
    //    irelcicr = grelcicr * (jsr - cai);
    //    tcicr = tcicr + dt;

    //    greljsrol = grelbarjsrol * (1.0 - Math.exp(-tjsrol / tauon))
    //            * Math.exp(-tjsrol / tauoff);
    //    ireljsrol = greljsrol * (jsr - cai);
    //    tjsrol = tjsrol + dt;

 
    //    djsr = dt * (itr - irelcicr - ireljsrol);
    //    bjsr = csqnbar - csqn - djsr - jsr + kmcsqn;
    //    cjsr = kmcsqn * (csqn + djsr + jsr);

    //    jsr = (Math.sqrt(bjsr * bjsr + 4.0 * cjsr) - bjsr) / 2.0;
        bjsr = 1 / (1 + (csqnbar * kmcsqn) / (Math.pow((kmcsqn + jsr), 2)));
        djsr = dt * bjsr * (itr - irel);
        jsr = djsr + jsr;
    }

    /**
     * Calculates new subspace ion concentrations
     */
    protected void concentrationSs() {
    bcass = 1 / (1 + (bsrbar * kmbsr) / (Math.pow((kmbsr + cass), 2)) + (bslbar * kmbsl) / (Math.pow((kmbsl + cass), 2)));
    idiffca = (cass - cai) / taudiff;
    dcass = -dt * bcass * ((ilca - 2 * inacass) * (acap / (vss * zca * frdy)) - irel * (vjsr / vss) + idiffca);
    cass = dcass + cass;
    //if (cass <= 0) cass = 0;

    idiffna = (nass - nai) / taudiff;
    dnass = -dt * ((ilcana + 3 * inacass) * (acap / (vss * frdy)) + idiffna);
    nass = dnass + nass;
    //if (nass <= 0) nass = 0;

    idiffk = (kss - ki) / taudiff;
    dkss = -dt * (ilcak * (acap / (vss * frdy)) + idiffk);
    kss = dkss + kss;
    //if (kss <= 0) kss = 0;
    }
}
    /**
     * Calculates new myoplasmic Ca ion concentration
     */
    //protected void concentrationCai() {
    //    dcai = -dt
    //            * (((caiont * acap) / (vmyo * zca * frdy))
    //                    + ((iup - ileak) * vnsr / vmyo)
    //                    - (irelcicr * vjsr / vmyo) - (ireljsrol * vjsr / vmyo));
    //    trpn = trpnbar * (cai / (cai + kmtrpn));
    //    cmdn = cmdnbar * (cai / (cai + kmcmdn));

    //    catotal = trpn + cmdn + dcai + cai;
    //    bmyo = cmdnbar + trpnbar - catotal + kmtrpn + kmcmdn;
    //    cmyo = (kmcmdn * kmtrpn) - (catotal * (kmtrpn + kmcmdn))
    //            + (trpnbar * kmcmdn) + (cmdnbar * kmtrpn);
    //    dmyo = -kmtrpn * kmcmdn * catotal;

    //    gpig = Math.sqrt(bmyo * bmyo - 3.0 * cmyo);

    //    cai = (2.0 * gpig / 3.0)
    //            * Math
    //                    .cos(Math.acos((9.0 * bmyo * cmyo - 2.0 * bmyo * bmyo
    //                            * bmyo - 27.0 * dmyo)
    //                            / (2.0 * Math.pow((bmyo * bmyo - 3.0 * cmyo),
    //                                    1.5))) / 3.0) - (bmyo / 3.0);
    //}
